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Abstract 

We propose a new method for the numerical solution of a PDE-driven model for colour image 
segmentation and give numerical examples of the results. The method combines the vector-valued 
Allen-Cahn phase field equation with initial data fitting terms. This method is known to be 
closely related to the Mumford-Shah problem and the level set segmentation by Chan and Vese. 
Our numerical solution is performed using a multigrid splitting of a finite element space, thereby 
producing an efficient and robust method for the segmentation of large images. 

This work has been submitted to the IEEE for possible publication. Copyright may be trans- 
ferred without notice, after which this version may no longer be accessible. 

I. Introduction 

The Mumford-Shah functional was first proposed in [24] as a general way to pose the problem 
of image segmentation. Reviews can be found for example in Petitot [26], and Fusco [12]. The 
generality of its statement has led to several methods of solution; in particular, the model sometimes 
referred to as the reduced Mumford-Shah was solved by the level set method by Chan and Vese [8], 
based on a previous paper on the motion of multiphase junctions tracked by the level set method 
by Zhao, Chan, Merriman and Osher [33], and subsequently extended in Chan and Vese [6], [7], 
[9], Vese [29], Chan, Shen and Vese [5]. Due to the extent of their work, it is also often referred 
to as the Chan- Vese model. The level set is used to track the boundaries of objects and should 
converge to the set of contours in the image. 

Esedoglu and Tsai [11] proposed the Allen-Cahn equation, also known as the phase field model, 
as a method of solution to the reduced Mumford-Shah problem when used in conjunction with the 
Chan- Vese fitting terms. The Allen-Cahn equation has been used in the context of image processing 
by Benes, Chalupecky and Mikula [3], who first convolve the image with a Gaussian smoothing 
kernel to eliminate noise, and then use it as an anisotropic gradient filter based on the proposals 
by Perona and Malik [20]. Our approach differs significantly in that we use no smoothing kernel 
and we modify the energy functional by adding fitting terms; moreover, we use the vector-valued 
formulation of the Allen-Cahn equation due to Garcke, Nestler and Stoth [13]. 

A different phase transition model (Modica-Mortola) was also recently used by Jung, Kang and 
Shen [14] , and although the model and numerical method of solution therein differs from this work, 
it is in many ways similar in the fundamental approach of adapting a phase transition model to 
solve the Mumford-Shah problem, and in the results obtained. 
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Our computations are solved by a multigrid algorithm which falls into the category of Successive 
Subspace Corrections (see Xu [30], [31]). This was successfully applied to the vector- valued Allen- 
Cahn equation in Kornhuber, [17], Kornhuber and Krause [18], [19] with a small variation (see 
Kornhuber and Krause [15]). 

In section [II] we briefly introduce and summarise previous directly relevant work leading up to 
section III-C] in which we formally introduce our own formulation and show how the minimisation 
of our functional leads to the desired system of PDEs; in section ITO1 we discretise the system and 
introduce the numerical method of solution, and in section IIVI we present a few practical aspects 
of implementation together with examples. 



A. Relation to the Mumford- Shah functional 

Given an image f 6 Q c I 2 , the Mumford-Shah method seeks to partition the domain Q into 
several subdomains Q; separated by a set K of boundaries, also known as edges or discontinuities. 
This segmentation takes the form of piecewise smooth functions u € Q — K which are discontinuous 
on K; the method of selection from all possible functions u is minimisation of an energy functional 



where \i, A are positive constants. 

The first term minimises the variation of u and promotes its smoothness, the second term 
minimises the length of interfaces and determines the boundaries between Q,, and the third term, 
sometimes referred to as the fidelity or fitting term, minimises the variation between u and I. 

As noted in Petitot [26] , the coefficients fj. and A define several scales of the problem: low /i leads 
to fine-grain segmentation, high p. to coarse-grain results. Sensitivity to contrast is measured by 
(4A 2 /,() 1 ^ 4 and robustness to noise depends on A/i. 

Many variations on this theme have been proposed since its first formulation. Mumford and Shah 
themselves pointed out that a reduced form of the problem, referred to as the minimal partition 
problem, is the restriction of u to piecewise constant functions, i.e. u = Ci with each q a constant 
on each connected region Q,. The minimising values are then clearly the averages of J across each 
region. 

In the level set case, using by way of example a single function (p to segment an image containing 
only one object against a background, the Chan-Vese functional introduced in [8] replaces the 
fidelity term in MS(u) by two fitting terms, 



where C\, Ci are the average of I inside and outside (p = 0, respectively. Considering for a moment the 
ideal situation in which the image contains only one object, i.e. is split into two regions of roughly 
constant value clearly separated by a gradient boundary, these two terms are clearly minimised 
when the set (f) = coincides with the contour of the object, i.e. the set K. 

B. A phase-field formulation 

The Allen-Calm PDE was introduced in [1] to model the domain coarsening occurring after a 
phase transition. It follows the evolution of a function u(x) known as the order parameter, which 
smoothly varies between the values of and 1 across an interfac^ to represent which parts of the 

lr The original model was defined on the interval [-1,1], but it is convenient to consider [0,1] for our purposes, 
without loss of generality. 



II. Image segmentation by the Allen-Cahn equation 




(1) 
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material are in one phase or another. It is obtained by minimising the following energy functional: 

AC(u) = J e\Vu\ 2 + j£U(u) dx (3) 

The function ^¥(11) represents a potential that attains minimal values at the two extreme values of 
u. This is not in general a convex function, nor is it necessarily smooth. Esedoglu and Tsai [11], 
for example, chose the quartic double- well form of the potential, ^(w) = u 2 (l — u) 2 . 

Comparing the first two terms in ([1]) and ([3|), the first term is identical up to a scaling constant, 
while the role of the second term is quite similar since the parameter e is directly related to the 
width of interfaces between phases; it is well known that minimising W(u) as above reduces both 
the width and length of boundaries. In this paper, we examine the results of extending ([3]) to its 
vector- valued formulation and combining it with fitting terms such as those in ([2]). 

It is reasonable to suggest that the results obtained by a level set method and a phase-field 
method should be closely comparable because both are known to be equivalent to curve motion by 
mean curvature; for an overview, see for example [10]. 

The method of solution described in [11] follows the MBO thresholding scheme by Merriman, 
Bence and Osher [21], [22], which assigns to the order parameter either one or the other extremal 
value at each step; we propose to use the formulation known as the double obstacle instead, in 
which, W takes the form 

W(u) = O(u) + Q(u), 

where 

O(m) = 

is known as the indicator function on the set [0,1], and Q(u) is the concave quadratic 

Q(u) = u(l - u). 



u € [0, 1] 
+00 u £ [0, 1] 



C. A modified vector-valued Allen- Cahn equation 

Our primary objective is to achieve a fast, robust image segmentation method. Given a domain 
QcK 2 and an image I: Q — > W with c data channels or colours, we follow the motion of a function 
u with several (N) components that we want to adapt to the significant features of I. We propose 
to do this by solving a system of PDEs on a finite-element space by a multrigrid method, and we 
derive this system by minimising an energy functional of the form 

6 = f £ El^ + u ^ ~ m) + o(u) + Au ■ F(c, I) dx, (4) 
Jq 2 6 



where 



F(c,I) = (I - c) 2 , c = ^— , (5) 



1 



the quantity c representing the average of I in u, in other words being a measure of the oscillation 
of the data over the support of u. 

In order to achieve a simultaneous segmentation of I into arbitrarily many pieces, we refer to the 
vector-valued formulation of the Allen-Calm system was introduced in Garcke, Nestler and Stoth 
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[13], which allows one to consider an arbitrary number of components to the order parameter, now 
described by a single vector-valued function u e r V N in the function set defined as 



V N := Iv e (H 2 (Q)) N : vfc) = 1 a.e. in Q, v { > a.e. in Q 1 , 



(6) 



In other words, the vector-valued function u must pointwise remain on the N-dimensional Gibbs 
Simplex 



G N : = he 



N 



< xt 



(7) 



which is itself an (N — l)-dimensional subset of the hyperplane 



-N. . 



[x e 



»N 



N 

E- 

i=i 



1 



Trivially, x ; < 1 Vx € G N (though not so for all x 6 E N ). 

An N-dimensional extension to the potential function W(m) is now required, with N minima given 
by Uj = \,Ui±j = 0; the form 



N 



Q(u) = Ui(l - m) 



i=i 

has been used in the following, with O(w) now the indicator function on G N and W(u) = O(u) + Q(u) 
as for the standard double-obstacle Allen-Cahn. 

Following standard procedure, we wish to minimise the energy functional dH) to derive a pde to 
which we can introduce a pseudo-time stepping to find a global minimum. To find the variational 
derivative of ([!]) in a direction v, care must be taken to remain in G N ; in other words, we do not 
want our function to leave the allowed set. If considering S(u + av) for some a € M, it must be that 
Yti u + av = Yji u = 1, and hence JijCcv = 0. To that end, Garcke, Nestler and Stoth [13] introduce 
the hyperplane 



TL N : : 



'ue 



which is at all points tangent to E (and hence to 
defined by 



Yj Ui 

i'=l 







together with the projection operator T 



Tx: = x- —(x ■ 1 N ) • ljv 



(8) 



as acting on a vector x € M N , where In: = (1, 1 ... 1) £ M N . Geometrically, the two hyperplanes E 
and TE N are parallel; TE N passes through the origin; the vector 



N 



iN 



is normal to TE N and represents the shortest distance between E N and TE N . As N grows larger, 
this distance grows smaller. By construction, u — v £ TE N Vw, z> e E N . 

We now turn to the question of minimising (j3J). As is usual for the Allen-Cahn functional ([3]), we 
use a gradient descent method, i.e. we find the directional derivative of and obtain a variational 
inequality to be solved numerically. Using the notation (•, •) to indicate the standard inner product, 
and using the the N-subgradient 



dO(u): = k e M N | ®(v) - 0(u) > Z{v - u)\ , 
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since we have 

<£, u - v) > Vn,»e Jom O, V£ € dO(u). (9) 

it is easy to see that a simple minimisation with respect to the order parameter u, which subse- 
quently determines the constants C; as well (as seen in [8]), leads to 

~\ i 1 

<— S(u),Tv) > (—eAu — —u + AF(c,I) -I — £,,Tv) 
du e e 



and hence, by introducing a pseudo-time parametrisation, we have the inclusion 

(u t -eAu + T^~u + AF(c,I)j,v-u) 3 (-d<$(u),v - u) (10) 

> (11) 

Vu, V € < V N . An approximation to this variational inequality can naturally be sought in terms of a 
finite element method, as described below. 

III. Discretisation and Numerical Solution 

It has been shown that iterative solvers such as the Jacobi, Gauss-Seidel, Successive Over- 
Relaxation and multigrid methods can be reduced to a class known as Subspace Correction methods; 
these first appeared in Xu [30], see also Xu [31], Kornhuber [16]. Convergence of subspace correction 
methods has been examined for example in Tai and Xu [27] for convex optimisation problems and 
in Neuss [25]. 

In essence, it is shown that a viable alternative to constantly projecting the solution from one 
level to another is to once and for all project all basis functions from all multigrid levels onto 
the coarsest one, solving the problem by iterating through all basis functions. An example of this 
projection method is shown in figure [21 This method was successfully applied to the vector- valued 
Allen-Cahn equation by Kornhuber [17], Kornhuber and Krause [18], [19]. 

In subsection IIII-AI we introduce the finite element method to establish the notation; building 
upon that basis, subsection IIII-BI shows the multigrid discretisation used, while subsection IIII-CI 
details a few aspects of the numerical implementation due to the Gibbs Space constraint and the 
double-obstacle method. 



A. Finite Element Notation 

The continuous domain Q is split into a set of subdomains T , referred to as a triangulation, 
given by the set of triangles t such that 

teT 

The natural length scale associated with each triangulation is 

h: = max diam{x). 

teT 

The vertices of all triangles form a set of n points or nodes, and for the purposes of this work, each 
triangulation is further required not to have any hanging nodes, i.e. nodes that are not corners of 
a triangle. Each node j:, 6 Q, i = 1 . . . n is assigned a function t]i{x) such that 




1 if x = Xi) 
if x + Xi. 
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Although these functions can be arbitrarily elaborate while still satisfying the specified require- 
ments, continuous piecewise linear functions are more than adequate for second-order problems 
such as ours. The set of all such functions, which can also be written as 

S: = [r] € C(Q): T]\ T is linear Vt e T) , (12) 

forms a basis for the finite element space 

<V* = span \m\ n i=l (13) 

such that all functions u 6 < V h can be represented as 

n 

u\x) = Y u u)m{x) (14) 
i=i 

and such that continuous functions u on the original domain can be represented on the triangulation 
by their piecewise linear interpolant 7i:C(Q) — > *V . 
We use the lumped mass and stiffness matrices 

n 

M. = Y j M = (i li/ l), A:= (V ni ,V nj ) Vr] (/ r] ; 6^. 

;'=i 

To discretise the inequality in (jlip by the finite element method, we pass to the the weak 
formulation 

(u\ + T [~^u h + AF(c, I)) , v h - u h ) - e{Vu h , V(v h - u h )) > Mv h e <V h (15) 

so that the discrete solution u h is only required to have H 1 regularity. We can then consider 
the above problem as a series of sub-problems over each basis function. This is known to be 
equivalent to an iterative solver, such as a Gauss-Seidel method; the precise sequence in which 
these are considered can be altered for instance to follow a red-black Gauss-Seidel pattern. By 
simple application of the properties of the projection operator T, using the discretisation method 
outlined above, we have n problems of the form 

u t + T (-^ M + AF(c,I)JjM(v - u) - euA(v - u) > 0. (16) 

with some appropriate time discretisation to follow. The inequality is due to the multi- valued nature 
of the subgradient at the boundaries of G N ; each iteration in the numerical method is performed 
as though ([TBI) were a strict equality, and if the result lies outside the acceptable space, then it is 
projected appropriately, as described in section IIII-C1 

B. Multigrid Solution 

The rate of convergence of Gauss-Seidel methods requires 0(n) iterations per time-step, and 
0(n 2 ) computations overall to converge to a solution, which is less than optimal. Multigrid methods 
are known to be more efficient, in most cases being of only 0(n) complexity overall; they rely on 
constructing several nested finite element spaces, usually by refining a coarse or macro triangulation 
T\ L — 1 times, eventually giving the fine triangulation 71, where 

Ti c 77+1. 

Figure Q] shows a simple example of one such refinement. Each triangulation is associated with 
its own finite element function space. Using projection and restriction operators, the solution is 
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Fig. 1. a) A coarse grid; b) the same grid, after one refinement. 



computed on one level, usually by a Gauss-Seidel solver or equivalent, and then passed to another 
level to be corrected. The theoretical basis of the multigrid method lies in showing that the error 
at each iteration can be considered to have several components of varying frequency, and that 
each time the problem is solved on a grid of natural spacing h, the error components that are of 
frequency h or higher are all reduced significantly, while those with a lower frequency are barely 
affected. The multigrid method significantly improves convergence by using several levels to deal 
with many components of the error in every iteration. 

A recent development in multigrid methods has been to consider iterative solvers as part of the 
Succesive Subspace Correction framework, which simply considers a minimisation over a sequence 
of function spaces such as (fT3j) , each defined as the set spanned by the basis functions defined on a 
different grid. In the context of finite elements, such a hierarchy of function sets is readily provided 
by the basis functions at each level of refinement. To be more precise, given the problem ,•: find 
u such that 

(u h t + t(-^u> 1 + AF(c,r)) f rii - u h ) - e(Vu h , V(r/; - u h )) > f| ; e <V N '> 

an SSC method applied to this case consists of simply solving all such 'Pi ; by looping through all 
nodes i and levels j and updating the solution at each step, in whatever appropriate sequence the 
specific method demands. 

Consider the finite element discretisation in (|16|) : we firstly discretise in time using a backward 
Euler scheme to obtain the fully discrete problem. At each time step, j, we use the SSC method 
to efficiently update the approximation u k as follows: starting at the coarsest level and moving to 
the finest in a standard multigrid w pattern, for every basis function rf L , i = \,1,...,n\ j on level L 
we update u k+1 = u + a' L where ai satisfies the inequality 

A (u k + a\ri - ui / 2 ■ \ 
M i + r(--„> + AF(cU)) 



- eA(u k + at]) > 



Once all basis functions and levels have been looped over, the iteration is complete; if the solution 
satisfies some prescribed error tolerance, it is accepted and becomes the new iterate u ] ; otherwise, 
the iteration is repeated using the computed iteration as the new starting point. 

All basis functions r\ of all multigrid levels are projected onto the finest grid. An example of this 
projection is given in figure [5J Thus each iteration is of the form 
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Fig. 2. a) A basis function denned on a coarse grid; b) the same basis function projected onto a finer grid 



M 



u k + arj - u* 



+ ^(-^w f + AF(c f ,/))J - eA(u k + arj) = 



5t 

before projection into the required function space. It is enough to multiply both sides by rf and 
re-arrange the above by simple algebra to obtain an equation for a: 

'2 



arj 1 (M - e5tA)i] = rj 1 M 



M fc + 6£T(-u f -AF(c f ,7)J +eAw /: J. (17) 



C. Gibbs Space Constraint Pseudocode 

In the vector-valued case, since N functions need to be updated simultaneously for each basis 
function, it is clear that a is also vector-valued and of size N. In order to preserve the constraint 

N 

1=1 

it must be that for each update 



N 

J> = 0. 

1=1 

This is actually taken care of by the choice of minimisation, i.e. by involving the projection operator 
T and only allowing search directions u + aTv. 

Secondly, for the functions to remain in the Gibbs space G N at all times, at each iteration of the 
numerical scheme we must ensure that the result still lies in G N , i.e. in addition to the restriction 
on the sum of a above, all components of u must be positive. This requirement may be less than 
straightforward to enforce on coarse grids, where several nodes are affected by the change in a; 
this is also mentioned in Kornhuber [17] and is due to the large support of the coarse grid basis 
functions after projection onto the finest grid; for that reason, a truncated multigrid algorithm may 
be more efficient than a full v-cycle. 

The pseudo-code for a general subspace correction is given below; it is assumed that for each 
basis function, a list of affected nodes has been drawn so that the least amount of work possible is 
being done. 

1) Evaluate a as per (|17p ; then, calculate a trial version of u. 

2) Check for each function w, whether any of its entries have become negative. If not, no further 
work is required. 
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3) Otherwise, we have a split of the functions W; into two sets, Q being the set of W; with negative 
(unacceptable) values, and P being the remainder. 

For each function in Q, we need to determine a new (Xj such that Uj remains between and 1, 
i.e. such that the lowest value is zero. This is easily done and we refer to this correction as j3. 
Since we know that J^a = must be satisfied, if we decrease a* then we need to increase all 
other ctj,] ± i by a corresponding amount. However, we can exclude those functions that are 
already negative. Therefore, if there are k functions in P and we decrease en, by an amount 
jS, then we need to decrease all by jS/fc. 
Do this for each function in Q. 

4) All functions in Q are now guaranteed to lie between and 1. However, the successive 
corrections f> applied to the functions in P may have caused some of those to become negative; 
hence, repeat from step 2, but only considering a restricted set of functions. Repeat this 
process until there are no more negative functions. 

Since at least one function is being removed from the set of all available functions at each 
repetition of this process, it is a procedure of maximum order N complexity. 



A. Determination of Natural Scales 

In principle, it is necessary to consider the possibility that the interface between two areas of 
interest be only one pixel wide; in other words, that two neighbouring pixels belong to different 
objects in the image. It is well known that the interface width of the double obstacle Allen-Calm 
equation is of 0(e) - see for example [10] and references therein. This immediately suggests a natural 
length scale for the problem: an e smaller than the size of a pixel in the given data would not make 
sense, and a larger e would only blur edges and corners. All our computations have been performed 
using this natural length scale. 

The value of e also suggests an estimate for a plausible A; if the value o = Ae were much less than 
1, the effects of the fidelity terms would be negligible; conversely, if it were too much larger than 1, 
the interface width would be reduced to such a slim margin that it could no longer be considered 
a diffuse interface. It was found in our computations that a value of roughly o e [10, 100] gave the 
most interesting results. 

However, it must be noted that the numerical discretisation of the Allen-Cahn equation can 
only give interface motion results of a satisfactory level of accuracy when there are several nodes 
to represent each interface, ideally 8 but in practice at least 4. In the context of a multigrid 
simulation, this can be achieved quite naturally by taking the given image's pixel structure as the 
reference grid and refining it an appropriate number of times - assuming uniform refinement to be 
the simplest case, where the number of nodes roughly doubles each time, this implies at least three 
refinements. It is also possible to coarsen the regular pixel grid a few times, depending on the size 
of the image. Thus the coarsest grid in the final setup will probably not correspond to the grid on 
which the data was defined. 

Another consequence of the need for a refinement process is the need to represent the image 
data on a finer grid than the one it was defined on. In the present model, the only link between 
the evolution of the function u to the image to be processed are the fidelity terms of the form 



IV. Implementation 




with the constants 




c = 
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approximating the average of I on suppu. The method of implementation of these terms is therefore 
critical to obtaining a good result. Moreover, a decision needs to be made as to what function space 
the image I is assumed to belong to. There are several candidates, such as C°°(Q) smooth functions, 
which make the mathematical analysis much easier and are generally derived by convolving the 
image data with a Gaussian; Lipschitz continuous functions, which include piecewise linear inter- 
polants; and functions of bounded variation. This distinction can be of great practical importance 
to a finite element method because depending on the chosen approach to node placement and 
fidelity term implementation it may be necessary to compare the values of u and J at points that 
do not correspond to given pixel values, or to points that lie exactly on the boundary between two 
or more pixels (e.g. corners). Which space the image belongs to ultimately depends on what kind 
of data one is examining, but in this case we have chosen to consider the image as a set of piecewise 
constant values over each pixel; this is done in order not to artificially blur any edges before any 
computations have even taken place. 

Together with mesh refinement, it is necessary to project the values of the original pixels to 
the fine grid. This is not done by interpolating the data values in any way, in order not to create 
spurious gradients. This is necessary because the fitting terms in our method rely on calculating 
the average value of I in each region; introducing new gradients also introduces small areas of 
different average value, which may be erroneously identified as new objects. Therefore, every newly 
created node is assigned exactly the same value as one of its neighbours. This preserves sharp 
discontinuities; it also leads to some staircasing of the data, but only on a relatively small scale 
that should be locked out by fixing e on the coarse grid. 

B. Data Projection 

If it is indeed necessary to project the values of the original pixels to the fine grid, this should 
not be by interpolating the data values in any way, in order not to create spurious gradients. This 
is necessary because the fitting terms rely on calculating the average value of I in each region; 
introducing new gradients also introduces small areas of different average value, which may be 
erroneously identified as new objects in their own right. Therefore, every newly created node is 
assigned exactly the same value as one of its neighbours, which preserves sharp discontinuities. 
In other words, the image data is taken to represent an underlying function I e BV of bounded 
variation. 

There are at least two distinct options for the implementation of these terms using a finite element 
model. The first is direct projection of the image data, node by node. The second is to note that 
the image function I is only weakly represented in the FE system, in the sense that only the value 
of its integral over a region is contained in the fidelity terms, and therefore in principle only the 
value of its integral on each element is required; projection can therefore be carried out triangle by 
triangle. This makes perfect sense as long as the triangulation is nested. In other words, the fidelity 
terms can be computed by either of the two following methods: either one uses the standard mass 
matrix, 

Mif.= (rn,Tij), 

„ \^ , , ( r M-I- u\ 2 
f = L m "['--tm-u) 

n=l x ' 
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Fig. 3. A clarification of two different data projection schemes 



or one computes the integral of I from first principles: 

N T„ 



=1 f=i 



AH- 



EM- 

Each method is associated with its own advantages, disadvantages, and computational costs. It 
is worth noting that the errors associated with each one decrease with each mesh refinement. The 
former can be thought of as projection by node and the latter as projection by simplex; examples 
are shown in figure O 



C. Post-processing 

Fig. 2£i shows an example of the post-processed solution. The green and blue functions indicate 
segmented regions as identified by the program; using these, the average of the data, osc I, is 
used in each to obtain the denoised data (red). This process is easily achieved by multiplying said 
average by the component itself, then summing all the resulting components; this is referred to as 
the composite. Recall that the measure of osc I on the support of each component U; is already 
known and used actively in the fitting terms driving the evolution (see ©). Further, because each 
component has values not identical to or 1, notably at each interface, it is useful to round all 
values to either extremum, in such a way that only one component is equal to 1 and all others 
are at any given point; in this way, segmented regions are defined more precisely. This naturally 
leads to the rounded composite, the advantage of which is shown in|3}D, a comparison of the errors 
given by the composite and rounded composite. 
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(a) 



(b) 




Fig. 4. a) A comparison of noisy data (black), segmented regions (green and blue), and resulting denoised signal 
(red); b) errors for the composite (blue) and rounded composite (green) segmentation. 




Fig. 5. a) noisy data (black) and rounded composite (red); b)original noise (blue) vs. recovered noise (green). 



Figured shows the noisy data (black) and the resulting rounded composite (red); for comparison, 
Fig. Ob shows the original noise (blue) and the recovered noise (green) obtained by subtracting the 
rounded composite from the initial data. 

D. Two Dimensional Examples 

Figure [6] shows the results of our method as applied to an image of concentric circles, with some 
noise. Figure [7] shows a more interesting case depicting what appears to be several overlapping 
geometric solids superimposed on a chequered background, for comparison with the results in 
Jung, Kang and Shen [14]. 

We also wish to examine the case of multi-channel data, the most obvious example of which 
is that of colour images. Although there is no reason why the proposed method could not be 
extended to images acquired at different wavelengths, for example by combining a visible night- 
time photograph with an infrared scan of the same area, the following arguments will concentrate 
on the case of colour images captured in the visible spectrum. 

It may seem natural to perform the processing on each channel individually. However, this does 
not take into account the fact that the information gathered from each channel is not disjoint but 
has a certain, possibly stronger, correlation across channels rather than within the space of each 
channel - in other words, the same pixel on each channel most likely corresponds to the same object 
being photographed, whereas there's no such guarantee for any two adjacent pixels in the same 
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Fig. 6. Segmentation of four concentric circles at [.25 .95 .55 .75] over a background at level .1, with added random 
noise of amplitude .05. The remainder shows a combination of error and noise, scaled to show all detail 



channelH Also, the problem remains of how to combine the resulting information, and this is not 
a straightforward course of action in the case of segmentation. 

Several methods to non-trivially process a colour image have previously been examined. The 
structure-texture decomposition method suggested by Meyer [23] was recently extended to colour 
images by Aujol and Kang [2] applying the G-norm to the RGB space. Another interesting approach 
is presented for example in Tang, Sapiro and Caselles [28], where the chromaticity values of each 
pixel are mapped onto the unit sphere and then processed using a diffusion-based filter; the authors 
consider both isotropic and anisotropic types. A rather similar approach seems to have been used in 
Yu and Bajaj [32], developed seemingly independently. It was argued in Kang and Shen [4] that a 
chromaticity-brightness (CB) filter outperforms a hue-saturation-value (HSV) filter in combination 
with a total variation method. 

We propose to use multi-channel information in our model by adapting the fitting terms in ([5]) 
to use the information from each of the colour channels If. 

, m k^l 

Ci,j{ui,lj) = — r 

Jo"' 



2 This depends on the field of interest. For example, in astronomy it may well be the case that a channel of 
information corresponds to a specific wavelength of light, in which case it is not only true that any one object is 
highly unlikely to appear across all wavelenghts, but the same absorption spectrum across different wavelenghts 
being captured at the same pixel most likely indicates two similar objects (e.g. gas clouds), one behind the other. 
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Fig. 7. Segmentation of a geometric composite 

3 

F(c irP Ij) = ^ \Ij - Ci^UiJj)^ dx 

This method can be followed regardless of how the multi-channel information is encoded, i.e. it can 
in principle be applied to any colour space, and any number of channels. For simplicity, the RGB 
space has been used in computations so far. Figures [8] and [9] were obtained using the RGB space. As 
can clearly be seen, the use of the RGB space cannot distinguish between a difference in colour and 
one in luminosity, potentially leading to somewhat unrealistic segmentations of complex real-world 
images; this could possibly be remedied by using the L*a*b* space and ignoring the luminosity 
component. 

V. Conclusion 

The method here presented is a flexible method for image segmentation and denoising based 
on the combination of previously consolidated work. It can either be fine-tuned to specific goals 
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or be allowed to run with almost no user input other than the initial image, which can define its 
own natural scale. Thanks to the multigrid formulation, computational costs will not spiral out of 
control. The method can also be further extended in several directions, including adaptive mesh 
refinement and adaptive time-stepping to further reduce computational time; it can deal with any 
desired colour space and any number of input and output channels, independent of each other. 

The Esedoglu-Tsai formulation involving MBO thresholding has the advantage of dealing well 
with the well-known unstable equilibrium value of 1/2; however, we would suggest that there are 
important dynamics motivated by the Chan-Vese fitting terms, especially at short time-scales, that 
are lost by the thresholding mechanism. Our multigrid solution seeks to retain the computational 
swiftness while capturing the full dynamics of the equations. 
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